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ABSTRACT 

We present temporal analysis of the three outbursts of the X-ray millisecond pulsar SAX J1808.4- 
3658 that occurred in 1998, 2000 and 2002. With a technique that uses the \ 2 obtained with an epoch 
folding search to discriminate between different possible orbital solutions, we find an unique solution 
valid over the whole five years period for which high temporal resolution data are available. We revise 
the estimate of the orbital period, P or b = 7249.1569(1) s and reduce the corresponding error by one 
order of magnitude with respect to that previously reported. Moreover we report the first constraint 
on the orbital period derivative, —6.6 x 10~ 12 < P or b < +0.8 x 10~ 12 s s _1 . These values allow us 
to produce, via a folding technique, pulse profiles at any given time. Analysis of these profiles shows 
that the pulse shape is clearly asymmetric in 2002, occasionally showing a secondary peak at about 
145° from the main pulse, which is quite different from the almost sinusoidal shape reported at the 
beginning of the 1998 outburst. 

Subject headings: stars: neutron — stars: magnetic fields — pulsars: general — pulsars: individual: 
SAX J1808.4-3658 — X-ray: binaries 



1. INTRODUCTION 

The X-ray transient SAX J1808.4-3658 was discovered 
in September 1996 when it exhibited an outburst de- 
tected by the BeppoSAX Wide Field Cameras (in 't Zand 
et al. 1998). The source showed type-I X-ray bursts 
that led to the identification of the compact object as 
a neutron star (hereafter NS) and to the derivation of 
a distance of about 2.5 kpc (in 't Zand et al. 2001). 
The source was found in outburst again in April 1998, 
when the high temporal resolution of the proportional 
counter array (PCA) on-board the Rossi X-ray Tim- 
ing Explorer (RXTE), allowed the discovery of coherent 
401 Hz pulsations, making this source the first known 
accretion-driven millisecond pulsar (Wijnands & van der 
Klis 1998). Timing analysis performed on data collected 
over the period 11 — 18 April 1998, allowed the deter- 
mination of the orbital parameters reported in Table |21 
(Chakrabarty & Morgan 1998, CM98 hereafter). SAX 
J1808. 4-3658 was again detected in outburst (and ob- 
served with RXTE/PCA) in January 2000 (van der Klis 
et al. 2000) and in October 2002 (see Wijnands 2004 for 
a review). 

Correct ephemerides are of capital importance for the 
study of this and similar systems, since these give us the 
possibility to monitor the evolution of the spin and the 
orbital period of the system. These, in turn, give impor- 
tant information on the mechanism of the acceleration 
of a NS to millisecond periods and on the evolutionary 
history of the system, respectively. Orbital ephemerides 
are also important for performing studies, for instance, in 
other wavelengths. In particular any search for radio pul- 
sations during quiescence (to test the hypothesis that the 
radio pulsar switches on during quiescence) needs cor- 
rect orbital parameters to perform coherent search for 
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pulsations and to improve by several orders of magni- 
tude the detection sensitivity. In this paper we apply a 
temporal analysis on the three most recent outbursts of 
SAX J1808. 4-3658 in order to improve the orbital pa- 
rameters of the system; we find a new orbital solution 
which appears to be valid from 1998 to 2002 (the entire 
period of the RXTE observations of SAX J1808.4-3658 
up to date). In §2 we present the timing technique we 
developed to obtain the revised estimate of the orbital 
period; in §3 we apply this technique to the available 
RXTE/PCA datasets of the 1998, 2000, and 2002 out- 
bursts of the source; in §4 we draw our conclusions. 

2. A METHOD OF TIMING FOR THE ORBITAL PERIOD 

Let us consider a binary system of "true" orbital pe- 
riod Porb, orbital separation A, and negligible eccen- 
tricity in which a source (hereafter the primary) or- 
bits the barycenter of the system at a radius Ai = 
Amil{m,\ + m2) (where mi = Mi/M© and — 
A^/Mq are the masses of the primary and its com- 
panion, the secondary, respectively). The orbital mo- 
tion causes delays in the times of arrival, i arr , of a sig- 
nal - emitted by the primary at the time t em - on the 
plane passing through the barycenter and perpendicu- 
lar to the line of sight: i arr = t cm + Z{t cm )/c, where 
Z{t) = Ai sin i sin [(2Tr/T> OTh )(t - T*)] is the distance of 
the primary from the plane defined above, i is the in- 
clination angle between the normal to the orbital plane 
and the line of sight, T* is the time of passage at the 
ascending node, and c is the speed of light. With an 
observational estimate of the orbital parameters a\ sini, 
P or b, and T*, we can correct the arrival times of photons 
emitted by the primary for the delays induced by the 
orbital motion. A first order correction gives: 
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The high temporal resolution RXTE data of SAX 
J1808. 4-3658 span a period of AXdata ~ 4.5 yr, from 
April 1998 to October 2002, during which SAX J1808.4- 
3658 performed A^ max ~ 19700 orbital cycles and n max ~ 
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57 billion spin cycles. In order to perform a timing anal- 
ysis over the whole ATdata we must unambiguously asso- 
ciate n, i.e. the number of elapsed spin cycles since the 
1998 outburst, to any given t m . Neglecting any error 
on the spin period estimate P sp in, induced by the errors 
in the orbital correction, this association would be possi- 
ble if the intrinsic error on the estimate of P sp in, crp sp i n , 

is SUCh that <7p spm ^ Pspin/^max ~ 4 X 1CT 14 s. CM98 

reported <7p S pin ~ 5 x 10~ 12 s, thus ruling out the possi- 
bility of making an overall timing analysis over ATdata- 
However, the large value of iV max suggested us to de- 
velop a method of timing of the orbital period P 01 - b m 
order to improve the orbital period estimate of CM98. 
Our idea is to estimate the time of passage at a given 
point in the orbit, e.g. the ascending node T* , for any 
given orbital cycle N, and to fit these times T*s vs the 
corresponding N to improve the estimate of P or b and to 
give a value (or an upper limit) for the orbital period 
derivative P or b- This timing can be done providing that: 

i) we can unambiguously associate N, i.e. the number 
of elapsed orbital cycles since the 1998 outburst, to any 
given T*. This is possible if the error on P or b, crp or b, is 
such that crporb < Porb/Amax ~ 0.4 s. CM98 reported 
fporb ~ 1 x 10~ 3 s, thus making this association possible; 

ii) we are able to estimate the T*s for any orbital cycle N 
with an error smaller than or* d = (cf-* +-^V 2 cp or b) 1//2 ~ 
Naporb obtained from the propagation of the CM98 er- 
rors in the "theoretical" computation of the T* rcd at a 
given N as T p * red (A) = T * + N x P orb , where T * is the 
time of passage at the ascending node at the beginning 
of the 1998 outburst, estimated by CM98. 

Our strategy to estimate the values of T*s relies on the 
fact that the correction of the photon arrival times with 
cq. Q induces a "smearing" in a pulse profile obtained 
by an epoch folding search performed on the corrected 
lightcurve because the corrections performed are affected 
by the errors in the estimate of the orbital parameters. 
Indeed, simple differentiation of the Doppler effect for- 
mula P spi „(i) = Pspin x (1 - v(t) / a)' 1 , where v(t) = Z(t) 
is the speed along the line of sight, demonstrates that the 
leading term in the error induced by eq. on P sp i n is 

^Pspin 

< | sin</>|P spin ai sini(27r/P or b) 2 -/V m axO'p orb , where 
<p = (27r/P or b)(£— T*) is the orbital phase. Since | sin0 < 
1 we have op spln < P spin ai sin i(27r/P orb ) 2 iV max c7p 01 . b . 
This means that folding a lightcurve of length ~ 1 hour 
~ P>rb/2 at a period close to P sp i n results into delays 
up to D max ~ 7roisini(27r/P or b)iV maac o-p orb . Adopting 
the CM98 estimates of the orbital parameters and the 
related errors we get P max ~ 3.2 x 10 -3 s ~ 1.3P sp i n for 
the 2002 data, which means that a significant smearing is 
is induced by <tt* ~ -/VmaxOporb, which represents our 
ignorance of the "true" time of ascending node passage 
in the 2002 data determined by the original error on P or b 
propagated through the factor N max . 

Therefore we can "experimentally" estimate the T*s 
in the 2000 and 2002 outbursts by correcting each ~ 
1 hour length observation with eq. in which several 
trial values of T* are adopted. We then choose the T* 
for which the corresponding corrected lightcurve gives a 
folded pulse profile with minimum smearing, obtaining 
an error in the determination of that particular T* that 



is expected to be smaller than the corresponding gt* 3 - 
Since a folded pulse profile obtained by folding a 
lightcurve corrected with the "true" time of ascending 
node passage (and therefore not smeared) has a bigger 
pulsed fraction than one obtained by folding a lightcurve 
corrected with a "wrong" time of ascending node passage 
(and therefore significantly smeared), the experimental 
estimate of each T* can be done by the following proce- 
dure: a) for a given observation we produce a number of 
corrected lightcurves applying eq. Q with different val- 
ues of the time of ascending node passage T t * rial ; b) we 
then perform on each lightcurve an epoch folding search 
for periodicities fitting the \ 2 v s spin period curve close 
to its maximum with a Gaussian profile, to derive the 
maximum x 2 value, Xmax; f° r the given value of T t * ial ; c) 
we plot these Xmax vs their corresponding T t * ial and we 
fit the part of this curve close to its maximum again with 
a Gaussian choosing, as our experimental estimate of the 
time of passage at the ascending node T*, the value of 
T t * ria | for which this Gaussian has a maximum. The un- 
certainties on T* are computed using post-fit residuals, 
i.e. the average fluctuation of the T* values vs. the or- 
bital cycle number TV with respect the average trend (see 
below). This is the more conservative way of estimat- 
ing the errors in the case of data with similar statistical 
significance as expected since we have used datasets of 
comparable length (see Tab. 1): indeed the post-fit er- 
rors are ~ 8 times bigger than the errors based on a fit 
of the top of the Gaussian. 

This technique supplies us with values of the epoch of 
ascending node passage Tjy that can be compared with 
the ones predicted using the original estimate of orbital 
parameters T* rcd (N) = T * + N x P orb . We therefore 
plot the delays AT^ = T^-T* ieA (N) versus N (we note 
that N is derived unambiguously as the integer part of 
(T^ - T *)/P orb under the assumption that \T^ - T * < 
P orb that we have also verified a posteriori) , fit the data 
with the relation 

AT£f =a + (3N + 1 N 2 (2) 

and finally obtain our improved orbital parameters as 

Tq = Pq* + a, <7t o = CT Q , Porb = Porb + /?, OP orb = <7/3) 
Porb = 2 7 /P rb, <J Pmb = (2/P rb)[(7/Porb) 2 ^ +CT 2 ] 1/2 - 

Although we devoleped this technique independently, we 
subsequentely found that a similar one had already been 
used by Paul et al. (2001) and Naik & Paul (2004). 

3. OBSERVATIONS AND ANALYSIS 

Throughout this paper we used the ToO public domain 
data of the PC A on-board RXTE (Bradt et al. 1993). 
In particular, we analyzed data from the PCA (Jahoda 
et al. 1996), which is composed by a set of five Xenon 
proportional counters operating in the 2-60 keV energy 

3 Note that, to estimate T* selecting the lightcurve with the 
minimum smearing, the <Tp sp i n (caused by the uncertainty in the 
knowledge of the true time of ascending node passage) must induce 
a genuine broadening of the folded pulse profile, and not delays 
more or less constant during each observation that can be compen- 
sated, in the folding search procedure, by the selection of a different 
spin period. Since the induced delays are not constant, but vary 
with time as sin <j>, the epoch folding search procedure must be per- 
formed over a time span that is a relevant fraction of the orbital 
period, so that sin</> varies significantly. This condition is verified 
by our sample of observations that last ~ 1 hour ~ P or b/2. 
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range with a total effective area of ~ 7000 cm 2 . We 
considered a sample (listed in Table Q of the available 
PCA observations of SAX J1808.4-3658 taken during the 
1998, 2000 and 2002 outbursts in which the spin modu- 
lation was particularly strong. For the temporal analysis 
we used event mode data with 64 energy channels and 
122/xs temporal resolution; the arrival times of all the 
events were reported to the solar system bary center. 

For the first 2002 observation (see Table ^) we pro- 
duced a set of lightcurves performing orbital correction 
adopting cq. (JIJ, with a\ sin? and P or b as determined by 
CM98. According to the method described above, we 
adopted different T* spanning a range around the pre- 
dicted value resulting from T* led (N G ) = T* + N G P olh , 
where N G is our initial guess of N: N G — INT{(to2 — 
T *)/P orb }, with < 2 the starting time of the first 2002 
observation. The T* range spans ~ 40 s in steps of 2 s 
as <7t* d ^ ^maxCporb ^ 20 s adopting the CM98 esti- 
mates. We then performed an epoch folding search on 
each corrected lightcurve, but, surprisingly, no pulsation 
was detected in any of them. 

Thus we extended the range of trial values for T* to 
explore the entire orbit until we obtained lightcurves for 
which the epoch folding search produced \ 2 peaks (typi- 
cally with maxima of ~ 3000 over an average value of ~ 
40) around periods close to P S pm- Following the method 
described above we then fitted the curve of the maxima 
of the x 2 (obtained from the epoch folding search on each 
lightcurve) vs the corresponding T* with a Gaussian and 
derived our best estimate of the time of ascending node 
passage for the first 2002 observation 2 . We then es- 
timated A02, i.e. the number of orbital cycles elapsed 
since T* as N 02 = INT{(T* q2 ~ T*)/P orb } = 19657 and 
therefore adopted a new estimate of the orbital period 
P Q n r e b w that is obtained from P n r c b w = (T^ 2 - T£)/N 02 . 
We checked the consistency of our provisional estimate 
of the orbital period P^™ demonstrating that iVo2 can 
be unambiguously associated to by verifying that 
INT{(T* Nq2 - T *)/P orb } = INT{(7% - T *)/P o n ™}. 

With this new estimate we repeated the procedure de- 
scribed above on each observation listed in Tabled with 
different T* spanning a range of 40 s around the new pre- 
dicted value resulting from T* CW (N G ) = T * + NqP™^, 
where N G = INT{(£-T *) /P or b}, with t the starting time 
of each observation. Encouragingly we always detected 
the pulsation, although with different strength (as de- 
duced from the different maximum x 2 ), by performing 
the epoch folding search on the corrected lightcurves. 
This allowed us to determine a set of 15 T^, one for 
each observation. 

In Fig. 1 (top and middle panel) we plot the delays be- 
tween the values and the values predicted adopting 
the CM98 estimates AT^ = Tfr - T* Icd {N) versus the 
number of orbital cycles elapsed since T *: where N — 
INT{(T* N - T *)/P orb }, and T p * red (A) = T * + TV x P orb . 
It should be noted that all the delays are < P orb which 
demonstrates a posteriori the consistency of our method, 
since we were able to associate unambiguously N to any 
Tff. Fitting these delays with eq. we find the im- 
proved orbital parameters and the first constraint on the 
orbital period derivative that are reported in Table El 

We finally corrected all the observations using eq. 
JQl with the improved value of the orbital period P or b 



and again we produced, for each observation, a set of 
lightcurves with different T* spanning a range around the 
predicted value resulting from T* led (N) = + P orb x N 
that allowed us to determine the new set of 15 T^. We 
finally computed the "post fit" residuals as 1Z(N) = 
T* N — T* red (iV) in order to check if the residuals were 
now compatible with zero. In Fig. 1 (bottom panel) we 
plotted 1Z(N) versus N to demonstrate the consistency 
of our method. 

4. DISCUSSION 

For the estimate of the T*s we have corrected the pho- 
tons arrival times with eq. lfT|l. adopting the estimate of 
P or b and x = aisini/c given by CM98. This seemed 
appropriate since we have demonstrated that, adopting 
the value of <7p orb and a x reported by CM98, the lead- 
ing term in o- Pspin depends only on o~ T* rcd ~ A max crp orb . 
However, at the end of our analysis, we found that our 
best-fit orbital period is incompatible with that reported 
in CM98, and this explains why the use of the ephemeris 
published in CM98 allows to recover the coherent pulsa- 
tions only during the 1998 outburst. We want to stress 
that a test of the correctness of our orbital solution is 
given by the fact that with our revised orbital parame- 
ters we can recover the pulsations in all the observations 
between 1998 and 2002. 

We also found a P or b 7^ (at ler confidence level) which 
means that in principle its effects on orbital parameters 
during ATdata ~ 4.5 yr must be discussed. Indeed the 
orbital period has varied and, because of the three out- 
burst that occurred during ATdata, m\ and 7712 varied 
as a consequence of mass transfer. Both these variations 
produce, because of third Kepler's law, a variation of the 
orbital separation along ATdata which reflects in a secu- 
lar variation of x. We therefore checked a posteriori that 
both the secular variation of x and P or b did not affect our 
determination of the orbital parameters of the system. 
By differentiating V rb(t) = V or b(t = *o) + Porb x(t-t ) 
and the third Kepler's law (assuming that the inclina- 
tion angle has remained constant) , we have demonstrated 
that the uncertainties on P or b and x in the 2002 outburst 
are always dominated by the uncertainties in their orig- 
inal determination as given by CM98. In particular we 
obtain that the uncertainty on x induced by mass trans- 
fer is of the order of 10~ 9 cZ| 5 m7 1 P6 times the uncer- 
tainty reported by CM98, as soon as the orbital period 
derivative is |P or b < 3 x 10~ 12 , the number of orbital 
cycles elapsed since the 1998 outburst is N < iV max , and 
the averaged mass loss rate by the secondary during the 
whole period ATdata is not much greater than the average 
mass transfer rate during the three outbursts, which is 
Moutb/Cltr^Mgyr- 1 ) = 3.7 d^ r,nq l R 6 (as can be eas- 
ily determined from an analysis of the X-ray luminosity 
during the 1998, 2000, and 2002 outbursts), where d 2 . 5 
is the source distance in units of 2.5 kpc and R§ is the 
NS radius in units of 10 6 cm. Therefore, our assumption 
of using a% sin % and P or b as given by CM98 in eq. JIJ is 
justified a posteriori, since their associated uncertainties, 
even taking into account the effects of orbital evolution, 
are always comparable with those estimated by CM98. 

The improved estimate of P or b derived in this pa- 
per (whose associated error, <7p orb , is 10 times smaller 
than that reported by CM98) allowed us to correct a 
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lightcurve using eq. Q with the improved value of the 
orbital period P or b and the predicted value of T* result- 
ing from T* red (iV) = T* + P orb x N. As discussed in the 
fourth paragraph of §2, folding this lightcurve modulo 
Pspin for any time interval ATfoid > P or b/2 introduces 
delays in the reconstruction of the pulse profile up to 
Anax ~ 2ai sini(27r/P orb )iV max crp orb (since | sin0| < 1 is 
a periodic function of period P or b/2 and average value 
2/tt ~ 0.6). Adopting our new estimate of <r_p orb we 
get Anax ~ 2.0 x 10 _4 s ~ 0.08P spin for the 2002°data. 
Which means that if the number of bins in the pulse pro- 
file is up to few tens, we do not expect a significant dis- 
tortion because of the uncertainties in the orbital param- 
eters. In Fig. Elwe compare two pulse profiles obtained 
by folding modulo P sp i n (the CM98 estimate) ~ 20 and 
~ 8 hours of data of the 1998 and 2002 outbursts, respec- 
tively, corrected with eq. It is evident that the 2002 
pulse profile has significantly varied with respect to the 
almost sinusoidal shape showed in the 1998 outburst; in 
particular we find a secondary peak at about 145° from 
the main pulse. This variability of the pulse profile sug- 
gests that strong caution should be taken in carrying out 
the standard pulse timing analysis when computing the 
delays of the time of arrival of a pulse by cross-correlating 
the locally folded pulse profile with a sinusoid, as done 
in CM98. Temporary variations of the pulse shape could 
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introduce spurious delays that may significantly affect 
the robustness of some derived parameters e.g. the local 
spin period derivative. 

In our analysis we found —6.6 x 10~ 12 < P or b < 
+0.8 x 10~ 12 s s" 1 (90% c.l.). Indeed, orbital evolu- 
tion calculations show that the magnitude of the P or b 
caused by angular momentum losses associated with the 
emission of gravitational radiation is expected to be one 
order of magnitude smaller, P or b ~ —3.0 x 10~ 13 (n — 

l/3)mi(mi + m2)~ 1 / 3 P 2h 2 ^ 3 s s -1 , where we assumed a 
mass-radius relation for the secondary R2 oc m^, and 
P2h is the orbital period in units of 2 hours. However 
the error in the estimate of P or b scales as AT ( | ata . This 
means that if, as expected, the next outburst of SAX 
J1808. 4-3658 will occur in 2005, we will reduce ap r by 
a factor ~ (4.5yr/7yr) 2 ~ 0.4. This new measure could 
severely constraint some of the evolutionary scenario pro- 
posed: for instance, the degenerate nature of the com- 
panion (that has been proposed to be a brown dwarf by 
Bildsten & Chakrabarty 2001) will be effectively probed, 
since n < 0, as expected for a degenerate companion, 
would imply P or b > 0. 
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TABLE 1 

RXTE OBSERVATIONS OF SAX J1808. 4-3658 USED in this work 



Observation ID Time of observation Phase 1 Phasc(CM) 2 



1998 outburst 



30411-01-01-02S 50914 19:26 - 19:52 0.18799(8) 0.18799(8) 
30411-01-03-00 50919 17:21 - 18:25 0.74694(8) 0.74724(8) 
30411-01-05-00 50921 07:48 - 08:57 0.84090(8) 0.84130(8) 



2000 outburst 



40035-05-03-00 51576 02:30 - 03:29 0.90879(8) 0.94999(8) 
40035-05-05-00 51582 00:36 - 01:29 0.47181(8) 0.51338(8) 



2002 outburst 
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Orbital phase of the beginning of each observation, referred to 
T* = 50914. 878465(l)Mjn and derived using our estimate of P orb = 
7249.1569(l)s 

2 Orbital phase of the beginning of each observation, referred to Tq — 
50914. 878465(1)MJD and derived using CM98 estimate of P orb = 
7249.119(l)s 



TABLE 2 

Improved orbital parameters of SAX J1808. 4-3658 



CM98 This work 



P orb (s) 7249.119(1) 7249.1569(1) 

A>rb (1CT 12 s s- 1 ) -6.6 < P orb < +0.8 

T * (MJD) a 50914.878465(1) 50914.878468(4) 



Note. — Numbers in parentheses are the 1 a uncertainties in the 
last significant figure. Limits on P or b are quoted at 90% confidence 
level. 

a CM98 reported a value of 50914.899440(1) MJD for the epoch of 
superior conjunction, i.e. when the NS is behind the companion; as, 
in this work, we considered the epoch of ascending node passage as 
a reference time, the CM98 reference time reported here has been 
decremented by P or b/4. 
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Fig. 1. — Upper panel: Delays, AT^, with respect to the ephemeris of CM98 vs the number of orbital cycles N elapsed since April 11, 
1998. Note that marks are much larger then error bars. Middle panel: expanded view of the three regions where points are clustered. 
Bottom panel: residuals 1Z(N) vs N. 
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Fig. 2. — 64-bin folded pulse profiles (P sp i n = 2.4939197575 ms, as reported by CM98). Upper panel a) refers to a folded light curve from 
18-04-1998 03:12:01 to 19-04-1998 00:57:43. Lower panel b) refers to a folded light curve from 17-10-2002 03:59:35 to 17-10-2002 11:59:35. 



